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Abstract 

The past two decades showed a rapid growing of 
physically-based modeling of fluids for computer graphics 
applications. In this area, a common top down approach is 
to model the fluid dynamics by Navier-Stokes equations and 
apply a numerical techniques such as Finite Differences or 
Finite Elements for the simulation. In this paper we focus 
on fluid modeling through Lattice Gas Cellular Automata 
( LGCA )for computer graphics applications. LGCA are dis- 
crete models based on point particles that move on a lattice, 
according to suitable and simple rules in order to mimic 
a fully molecular dynamics. By Chapman-Enskog expan- 
sion, a known multiscale technique in this area, it can be 
demonstrated that the Navier-Stokes model can be repro- 
duced by the LGCA technique. Thus, with LGCA we get 
a fluid model that does not require solution of complicated 
equations. Therefore, we combine the advantage of the low 
computational cost of LGCA and its ability to mimic the 
realistic fluid dynamics to develop a new animating frame- 
work for computer graphics applications. In this work, we 
discuss the theoretical elements of our proposal and show 
experimental results. 



1. Introduction 

Physically-based techniques for the animation of natu- 
ral elements like fluids (gas or liquids), elastic, plastic and 
melting objects, among others, have taken the attention of 
the computer graphics community 1 14|. The motivation for 
such interest rely in the potential applications of these meth- 
ods and in the complexity and beauty of the natural phe- 
nomena that are involved 1251 l3l. In particular, techniques 
in the field of Computational Fluid Dynamics (CFD) have 
been applied for fluid animation in applications such as vir- 
tual surgery simulators, computer games and visual effects 

(UEIl. 

In this paper we focus on physically-based fluid anima- 
tion for computer graphics applications (see [Q and ref- 



erences therein). Basically, the works in this area fall in 
to two categories: Realistic fluid and Interactive, or Real- 
Time, fluid animation. The former is more suitable for the 
special effects industry fvf\ while the later is appropriate 
for interactive applications like computer games and vir- 
tual surgery l23l 1151 . The work (Sj is a remarkable one in 
this area which includes fluid equations and numerical tech- 
nique |9|, shortly Computational Fluid Dynamics (CFD), 
and scientific visualization methods 1201 . The literature of 
this field reports gas f5','251 and water simulations 1 13 1, in- 
teraction between liquids and deformable solids 116J . and 
ofliers l24im. 

A majority of fluid animation methods in computer 
graphics use 2D/3D mesh based approaches that are math- 
ematically motivated by the Eulerian methods of Finite El- 
ement (FE) and Finite Difference (FD), in conjunction with 
Navier-Stokes equations of fluids [9J. These works are 
based on a top down viewpoint of the nature: the fluid is 
considered as a continuous system subjected to Newton's 
and conservation Laws as well as state equations connect- 
ing the macroscopic variables of pressure P, density p and 
temperature T. 

In this paper, we change the viewpoint to the bottom 
up model of the Lattice Gas Cellular Automata (LGCA) 
(6^1 . These are discrete models based on point particles that 
move on a lattice, according to suitable and simple rules in 
order to mimic a fully molecular dynamics. Particles can 
only move along the edges of the lattice and their interac- 
tions are based on simple collision rules. There is an exclu- 
sion principle that limits to one the number of particles that 
enter a given site (lattice node) in a given direction of mo- 
tion. Such framework needs low computational resources 
for both the memory allocation and the computation itself. 
Such models have been applied for scientific application in 
two-phase flows description (gas-liquid systems, for exam- 
ple), numerical simulation of bubble flows (10|, among oth- 
ers. Besides, Wolfram |26| has studied the computational 
and thermodynamics aspects of these models for fluid mod- 
eling. 

In this paper we focus on fluid modeling through Lat- 



tice Gas Cellular Automata (LGCA) for computer graphics 
applications. Specifically we take a special LCGA, intro- 
duced by Frisch, Hasslacher and Pomeau, known as FHP 
model, and show its capabilities for computer graphics ap- 
plications. By Chapman-Enskog expansion, a known multi- 
scale technique in this area, it can be demonstrated that the 
Navier-Stokes model can be reproduced by FHP technique. 
However, there is no need to solve Partial Differential Equa- 
tions (PDEs) to obtain a high level of description. There- 
fore, we combine the advantage of the low computational 
cost of LGCA and its ability to mimic the realistic fluid 
dynamics to develop a new animating framework for com- 
puter graphics applications. Up to our knowledge, there are 
no references using FHP for fluid animation in Computer 
Graphics. In this work, we discuss the theoretical elements 
of our proposal and present some experimental results. 

The paper is organized as follows. Section |2] offer a re- 
view of CFD for fluid animation. Section |3l describes the 
FHP model its multiscale analysis. The experimental re- 
sults are presented on section |4] Conclusions are given on 
Section|5] 

2. Navier-Stokes for Fluid Animation 

The majority fluid models in computer graphics follow 
the Eulerian formulation of fluid mechanics; that is, the fluid 
is considered as a continuous system subjected to Newton's 
and conservation Laws as well as state equations connecting 
the macroscopic variables that define the thermodynamic 
state of the fluid: pressure P, density p and temperature T. 

So, the mass conservation, also called continuity equa- 
tion, is given by ||9l: 



| + V.(p.)=0 



(1) 



The linear momentum conservation equation, also called 
Navier-Stokes, can be obtained by applying the third New- 
ton's Law to a volume element dV of fluid. It can be written 
as ID.: 



du ^ 
-+U-WU 



1. 



VP + F + /i V^M + -V (V-u) 



(2) 

where F is an external force field and n is the viscosity of 
the fluid. Besides, the equation V-u = must be added 
to model incompressible fluids. Thus, if we combine this 
equations with expression Q we obtain the Navier-Stokes 
equations for incompressible fluids (water, for example): 



du ^ „^ 



-VP + F- 



(3) 



(4) 



Also, we need an additional equation for the pressure 
field. This is a state equation which ties together all of the 
conservation equations for continuum fluid dynamics and 
must be chosen to model the appropriate fluid (i.e. com- 
pressible or incompressible). In the case of liquids, the pres- 
sure P is temperature insensitive and can be approximated 
by P = P (p). Morris in 1121 proposed an expression that 
have been used for fluid animation also 1 1 3l : 



P 



2 
C p 



(5) 



where c is the speed of sound in this fluid 1211 . 

Equations (|3|l-(|5|l need initial conditions 

(p {t = 0, X, y,z) ,u{t = 0, X, y, z)). Besides, in practice, 
fluid domain is a closed subset of the Euclidean space and 
thus the behavior of the fluid in the domain boundary - 
boundary conditions - must be explicitly given. For a fixed 
rigid surface S, one usual model is the no-sleep boundary 
condition that can be written as: 



0. 



(6) 



Also, numerical methods should be used to perform 
the computational simulation of the fluid because the fluid 
equations in general do not have analytical solution. Fi- 
nite Element (EE) and Finite Difference (ED) are known 
approaches in this field. Recently, the Lagrangian Method 
of Characteristics 1221 1231 and the meshfree methods of 
Smoothed Particle Hydrodynamics (SPH) |13| und Moving- 
Particle Semi-Implicit (MPS) 1 19| have been also applied. 

If the fluid is temperature sensitive, then an energy con- 
servation law should be applied. For example, in 1 5 1 authors 
develop a framework for hot turbulent gas animation. The 
model comprises equations (|3},@,(|6j as well as the follow- 
ing equation for temperature change and the buoyant force, 
respectively: 



BT 

— = AV^P - V • (Tu) , 



(7) 



(8) 



where A is the diffusion coefficient, Tq is a reference tem- 
perature and /3 is the coefficient of thermal expansion. The 
numerical method used in |5| is Finite Difference. This 
work can reproduce a hot gas behavior with some realism 
but has the limitation that the integration time step is con- 
strained to: 



At < 



(9) 



where h is the mesh resolution. Besides, the restriction of 
equation Q is not suitable for a compressible system like a 
gas. 



Henceforth, after that work, we can find works that: (a) 
Propose more stable models to achieve faster simulations; 
(b) Use truly meshfree Lagrangian methods; (c) Include 
realistic behaviors of truly incompressible flow simulation 
and interaction of fluids with deformable solids; (d) Use 
GPU capabilities in order to achieve faster simulations for 
interactive applications; (e) Generate special effects through 
fluid flows; among others 1 1 1. 

From the viewpoint of fluid models, all the cited works 
are top down approaches in the sense that the relationships 
of interest are between variables that capture the global 
properties of the system; that is, pressure, density and tem- 
perature. These relationships are expressed in ordinary or 
partial differential equations like (|3}- 

On the other hand, bottom up models start from a de- 
scription of local interactions. These models usually in- 
volve algorithmic descriptions of individuals, particles in 
the case of fluids. Analysis and computer simulation of bot- 
tom up models should produce, as emergent properties, the 
global relationships seen in the real world, without these 
being built into the model. Thus, there is no need to use a 
PDEs and numerical methods to obtain a high level of de- 
scription. 

For Computer Graphics applications, such approach is 
explored in for real-time simulation and animation of 
phenomena involving convection, reaction-diffusion, and 
boiling. An extension of cellular automata known as the 
coupled map lattice (CML) is used for simulation. CML 
represents the state of a dynamic system as continuous val- 
ues on a discrete lattice. In II II the lattice values are stored 
in a texture, and pixel-level programming are used to imple- 
ment simple next-state computations on lattice nodes and 
their neighbors. However, Navier-Stokes models are not 
considered and CML still uses continuous values for rep- 
resentations. That is also the case of Lattice Boltzmann 
models |8 1. In this paper we propose the application of an 
even more simples model, the FHP one, for fluid simula- 
tion. It will be demonstrated how Navier-Stoke models can 
be reproduced by this method. FHP is described in the next 
section. 

3. FHP and Navier-Stokes 

The FHP was introduced by Frisch, Hasslacher and 
Pomeau fl] in 1986 and is a model of a two-dimensional 
fluid. It can be seen as an abstraction, at a microscopic 
scale, of a fluid. The FHP model describes the motion of 
particles traveling in a discrete space and colliding with 
each other. The space is discretized in a hexagonal lattice. 

The microdynamics of FHP is given in terms of Boolean 
variables describing the occupation numbers at each site of 
the lattice and at each time step (i.e. the presence or the 
absence of a fluid particle). The FHP particles move in dis- 



crete time steps, with a velocity of constant modulus, point- 
ing along one of the six directions of the lattice. The dynam- 
ics is such that no more than one particle enters the same 
site at the same time with the same velocity. This restric- 
tion is the exclusion principle; it ensures that six Boolean 
variables at each lattice site are always enough to represent 
the microdynamics. 

In the absence of collisions, the particles would move 
in straight lines, along the direction specified by their ve- 
locity vector. The velocity modulus is such that, in a time 
step, each particle travels one lattice spacing and reaches a 
nearest-neighbor site. 

In order to conserve the number of particles and the mo- 
mentum during each interaction, only a few configurations 
lead to a non-trivial collision (i.e. a collision in which the 
directions of motion have changed). When exactly two par- 
ticles enter the same site with opposite velocities, both of 
them are deflected by 60 degrees so that the output of the 
collision is still a zero momentum configuration with two 
particles. When exactly three particles collide with an an- 
gle of 120 degrees between each other, they bounce back to 
where they come from (so that the momentum after the col- 
lision is zero, as it was before the collision). Both two- and 
three-body collisions are necessary to avoid extra conserva- 
tion laws. Several variants of the FHP model exist in the 
literature |4|, including some with rest particles like mod- 
els FHP-II and FHP-III. For all other configurations no col- 
lision occurs and the particles go through as if they were 
transparent to each other 

The full microdynamics of the FHP model can be ex- 
pressed by evolution equations for the occupation numbers 
defined as the number, Ui (r, t), of particle entering site f 
at time t with a velocity pointing along direction Ci, where 
i — 1,2, ... ,6 labels the six lattice directions. The numbers 
n.i can be or 1. 

We also define the time step as At and the lattice spacing 
as Ar. Thus, the six possible velocities Vi of the particles 
are related to their directions of motion by 

V^ = -^C.. (10) 

Without interactions between particles, the evolution equa- 
tions for the fii would be given by 

n^{f+ArC^,t + At)^n^{f,t) (11) 

which express that a particle entering site r with velocity 
along Ci will continue in a straight line so that, at next time 
step, it will enter site r + A^q with the same direction of 
motion. However, due to collisions, a particle can be re- 
moved from its original direction or another one can be de- 
flected into direction ci. 

For instance, if only Ui and 71^+3 are 1 at site r, a 
coUision occurs and the particle traveUng with velocity Vi 



will then move with either velocity Vi^i or Vi+i, where 
i = 1, 2, . . . , 6. The quantity 

Di = niUi+s (1 - Hi+i) (1 - (1 - f^^+4) (1 - nj+s) ■ 

(12) 

indicates, when Di = 1 that such a collision will take place. 
Therefore ni — Di is the number of particles left in direction 
Ci due to a two-particle collision along this direction. 




Figura 1 : The two-body collision in the FHP. 



Now, when n, = 0, a new particle can appear in direction 
Ci, as the result of a collision between ri^+i and 71^+4 or a 
collision between n^-i e 71^+2- It is convenient to introduce 
a random Boolean variable q (r, t), which decides whether 
the particles are deflected to the right (q = 1) or to the left 
(q = 0), when a two-body collision takes place. Therefore, 
the number of particle created in direction Ci is 

gA-i + (i-g)A+i. (13) 

Particles can also be deflected into (or removed from) di- 
rection Ci because of a three-body collision. The quantity 
which express the occurrence of a three-body collision with 
particles rii, 71^+2 and 77^+4 is 

Ti = nini+2ni+4 (1 - 7ii+i) (1 - 71^+3) (1 - 71^+5) (14) 

As before, the result of a three-body collision is to modify 
the number of particles in direction Ci as 

ni-Ti+ Ti+3, (15) 

Thus, according to our collision rules, the microdynamics 
of a LGCA is written as 

n^{r + ArC,,t + At) = n, {r,t) + {n {r,t)) (16) 

where 17; is called the collision term. 

For the FHP model, fti is defined so as to reproduce the 
collisions, that is 

= - A + gA-i + (1 - g) A+i - + T,+3. (17) 



Using the full expression for _D, and Ti, given by the Equa- 
tions (I12> -( I14> . we obtain, 

(18) 

= -7li72j+2"-i+4 (1 - (1 - ni+3) (1 - 71^+5) 

+ 7lj+l77j+377i+5 (1 - Ui) (1 - 7li+2) (1 - "1+4) 

- 71^77^+3 (1 - 77j+l) (1 - 77^+2) (1 - n.1+4) (1 - 77i+5) 

+ (1 - g) n^+i77i+4 (1 - 77j) (1 - ni+2) (1 - nj+3) 

+ (1 - g) (1 - 77,+5) 

+ g77j+2"i+5 (1 - rii) (1 - Ui+i) (1 - 7ii+3) (1 - 77^+4) . 

These equations are easy to code in a computer and yield a 
fast and exact implementation of the model 

Until now, we deal with microscopic quantities. How- 
ever, the physical quantities of interest are not so much the 
Boolean variables Ui but macroscopic quantities or average 
values, such as, for instance, the average density of parti- 
cles and the average velocity field at each point of the sys- 
tem. Theses quantities are defined from the ensemble av- 
erage Ni if, t) = {rii {r, t)) of the microscopic occupation 
variables. Note that, Ni {r, t) is also the probability of hav- 
ing a particle entering the site r, at time t, with velocity 

In general, a LGCA is characterized by the number z of 
lattice directions and the spatial dimensionality d. In our 
case d = 2 and z = 6. Following the usual definition of 
statitical mechanics, the local density of particles is the sum 
of the average number of particles traveling along, each di- 
rection Ci 

z 

p{r,t) = Y,N^{r,t). (19) 

i=0 

Similarly, the particle current, which is the density p times 
the velocity field li, is expressed by. 

z 

p [f, t) u {r, t) = Y.'^.N, (r-, t) . (20) 

Another quantity which will play an importante role in the 
up coming derivation is the momentum tensor 11 defined as 

z 

Ila/J = 'Y^ViaVipNi {f, t) (21) 

i=0 

where the greek indices a and (3 label the d spatial com- 
ponents of the vectors. The quantity 11 represents the flux 
of the a— component of momentum transported along the 
/3— axis. This term will contain the pressure contribution 
and the effects of viscosity. 



The starting point to obtain the macroscopic behavior of 
the CA fluid is to derive an equation for the N- s. Averaging 
the microdynamics ( I16> yields 

N, (r + A^c,;, t + At)- N, (f, t) = {n, {n (f, t))) (22) 

where Jl; is the colHsion term of the LGCA, under study. It 
is important to notice that 57^ (n) has some generic proper- 
ties, namely 



E 



n,^Q e ^u^n, = 



(23) 



expressing the fact that particle number and momentum are 
conserved during the collision process (the incoming sum 
of mass or momentum equals the outgoing sum). 

The Ni's vary between and 1 and, at a scale L >> Ar 
e T >> Af, one can expect them to be smooth functions 
of the space and time coordinates. Therefore, equation (I22> 
can be Taylor expanded up to second order and gives 



A^ (c, • V) N, {r, t) + AtdtN, (r, t) 



(24) 



+ i (A,)' {c^ ■ V f N, (r, t) + ArAt (c, • V) dtN, (f, t) 

+ \{^tf idtfN,{f,t) = {n, inir,t))). 

where (9*)^ is the second derivative in respect to the time 
parameter t. 

At a macroscopic scale L >> Ar, following the proce- 
dure of the so-called multiscale expansion L18J . we intro- 
duce a new space variable fi such that 



ri — e9f J e dr = edp-^ 



(25) 



with e << 1. We also introduce the extra time variables 
ti — et and t2 — e^t, as well as new functions depend- 
ing on fi, ti and t2, Nf = (ti, ^2, ^i) and substitute into 
equation (I24> 



N,, m 



dt -> ec?ti + e^t^t2 9r (26) 



together with the corresponding expressions for the second 
order derivatives. Then obtain new equations for the new 
functions N^. Thus, we may write 1 18 1, 



m 



(27) 



The Chapman-Enskog method is the standard procedure 
used in statistical mechanics to solve an equation like (I24> 
with a perturbation parameter e. Assuming that (57^ (n)) 
can be factorized into 51^ (N), we write the contributions of 
each order in e. According to multiscale expansion (??, the 
right-hand side of ( I24> reads 



da 



(1) 



CI 



(28) 



Using expressions (I25> -( I27> in the left-hand side of ( I24t and 
comparing the terms of the same order in e in the equation 
(HD, yields 



and 



r(0) 



A,2-.[ ON, ) 



(29) 



(30) 



N, 



(1) 



where the subscript 1 in spatial derivatives (e.g. dia) indi- 
cates a differential operator expressed in the variable ri and 



^ {ci ■ Vri ) = diaVia, from equation ( fTol i. 



We also impose the extra conditions that the macroscopic 
quantities p and pu are entirely given by the zero order of 
expansion Ml\ 



(0) 



and 



pu 



(0) 



(31) 



and therefore 



^iVj'' = and J^^^^i " for ? > 1 (32) 

Thus, following the Chapman-Enskog method we can 
obtain (5), from equation J24> . the following result at order 
e 

dt^ p + divi pu = (33) 

and 



t(0) 







(34) 



dtiPUa + dll3^al3 

On the other hand, if we consider the terms of order and 
using the relations ( I33> and ( I34> to simplify, we have 



dt2PUa + d 



'1/3 



= 
(35) 

The last equation contains the dissipative contributions to 
the Euler equation i34l . The first contribution is li^^^ which 
is the dissipative part of the momentum tensor The sec- 
ond part, namely ^ (dti^al ~^ '^q'/37) comes from the 
second order terms of the Taylor expansion of the discrete 
Boltzmann equation. These terms account for the discrete- 
ness of the lattice and have no counterpart in standard hy- 
drodynamics. As we shall see, they will lead to the so-called 
lattice viscosity. The order e e can be grouped together to 
give the general equations governing our system. Summing 
equations i33i and i35\ with the appropriate power of e as 
factor and we obtain the continuity equation (see expression 



dtp + div pu ^ 



(36) 



Similarly, equation ( I34t and ( I35t yields (5] 



We can see ||6| that the lattice viscosity is given by 



dtpUa 



n, 



Q/3 



At 
2 



eduU 



(0) 



5 



(0) 



dr. "'^'^ 



(37) 

We now turn to the problem of solving equation ( I29> to- 
gether with conditions i3l\ in order to find A^^"' as func- 
tions of p and pil. The solutions A'^j*^^'' which make the col- 
lision term fl vanish are known as the local equilibrium so- 
lutions. Physically, they correspond to a situation where the 
rate of each type of collision equilibrates. Since the colli- 
sion time Af is much smaller than the macroscopic obser- 
vation time, it is reasonable to expect, in first approximation 
that an equilibrium is reached locally. 

Provided that the collision behaves reasonably, it is 
found [6j that the generic solution is 



(0) 



1 



1 + cxp (^-A - B - v^ 



(38) 



This expression has the form of a Fermi-Dirac distribution. 
This is a consequence of the exclusion principle we have 
imposed in the cellular automata rule (no more than one 
particle per site and direction). This form is explicitly ob- 
tained for the FHP model by assuming that the rate of direct 
and inverse collisions are equal. The quantities A e B are 
functions of the density p and the velocity field u and are to 
be determined according to equations i'iH . In order to carry 
out this calculation, iV^^"^ is Taylor expanded up to second 
order in the velocity field u. One obtains Q 

Nl ' ^ ap+ -irVi ■ u H j—QiapUaUp (39) 

where a, /3, 7 are summed over the spacial coordinates, e.g. 
a,/3,7e{l,...,d},i;= A^,a= i,6 = f and 



}ial3 — ViaVip 



(40) 



The function G is obtained from the fact that iv|°^ is the 
Taylor expansion of a Fermi-Dirac distribution. For FHP, it 
is found 12,,6J 

2 (3 - p) 



G{p) 



3 (6 - p) 



We may now compute the local equilibrium part of the 
momentum tensor, Tl^^l and then obtain the pressure term 



p = aC2V p - 



C2 



pG{p)u' 



(41) 



where C2 



^lattice — —Cib 



d{d + 2) z 2 



-A. 



2 (d + 2) 



t 2 

V 



The usual contribution to viscosity is due to the collision 
between the fluid particles is given by (6| 



Ucoll ^ AtV 



A 



where — A is given by —A = 2s (1 — s) where s = f 
Therefore, the Navier-Stokes equation reads 

dtu + 2C4G{p){u-V)u = --\/p + iyS/^u (42) 

P 



where 



= A.v'bG, ±-1]=^ A^l] (43) 



) - ( 






J d + 2^ 






our discrete 


fluid. 





4. Experimental Results 

In this section we describe some experiments with FHP 
for bidimensional fluid simulation. Firstly, we highlight the 
simplicity of creating new configurations. Figure 2 shows 
an initial configuration with zero density in the middle of 
the system. It is not required any extra mathematical ma- 
chinery to deal with such density discontinuity because sys- 
tem rules do not undergo modifications. Figures 1 were gen- 
erated with 80.000 particles with position and velocity di- 
rections randomly distributed. The lattice resolution is 100 
by 100 points. 

The density distribution at time 10 and 25 (Figures 3 and 
4) show an interesting pattern near the front of the discon- 
tinuity. The evolution for time 50 is even more interesting 
(Figure 5). If we want to predict such effects, we need to 
consider Navier-Stokes equations. However, if the aim is to 
explore the visual effect, we can just simulate and take the 
desired result at its time. As expected, the system evolves 
towards a configuration in thermodynamic equilibrium (or 
maximum entropy |26|). Figure 6 shows such state. From 
the macroscopic viewpoint, the fluid achieves a static con- 
figuration in which the macroscopic velocity u is null every- 
where. If we decrease particle density, the pattern obtained 
is basically the same, as we can verify through Figure 7. 

The configuration of pictured on Figure 2 can be gen- 
eralized by an initial density with a disconnected zero set. 
Figure 10-a pictures such example. We get an interesting 
pattern formation presented on Figure 10-b. These patterns 
evolve to the "S" formations pictured on Figure 1 1 . 



Besides, we can take advantage of the simplicity of the 
model for changing boundary. For a LGCA, there is no 
need to re-build the lattice. It is just a matter of finding the 
boundary cells of the lattice and apply the proper collision 
rules for particles entering the corresponding sites. Next, 
we show the tests using a homogeneous particles distribu- 
tion with velocity in the horizontal direction. It is interesting 
to observe the patterns at the right hand side of the Figure 
9. Particles that coUide with the domain boundary also will 
collide with the insident particles which increases the den- 
sity nearby. 

5. Conclusions 

In this paper we propose the FHP model for fluid model- 
ing in computer graphics applications. We discuss the the- 
oretical elements of our proposal and discuss some experi- 
mental results. Further works are the incorporation of ex- 
ternal forces and model two-fase systems for visual effects 
generation. 
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Figure 10: (a) Initial configuration, (b) Trasient pattern 
formation. 




Figure 1 1 : Evolution of the configuration pictured on Figure 9-a. 



